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INTRODUCTION 



A. BACKGROUND 

1. The Importance of Stealth 

The concept of stealth has been around since long before man made his appearance 
on earth. Nature realized that the survival of its creatures depended on their blending in 
with the background. Evolution has provided the means for nature’s creatures to develop 
and refine stealth capabilities. It has been during the last fifty or so years that the military 
forces of the world have taken a keen interest in stealth [1], The results were displayed 
by the success of the FI 17A stealth aircraft during the Persian Gulf War [2, 3, 4], 

The success of these aircraft during the Gulf War can be attributed to many 
factors. By far the most important factor was possession and use of stealth technology. 
The U.S. military had it and the enemy did not. More importantly, the enemy did not 
possess the technology to defeat a stealthy airborne platform. That was in 1991, more 
than six years ago. Since then stealth technology has been ever so slowly let out of its 
“black box.” It is no longer “our own little secret weapon” that only the U.S. military 
possesses, and there is no guarantee that during the next conflict our adversaries will not 
possess anti-stealth technology. In fact ultra-wideband radar (UWB) was being 
developed for detection of stealth platforms before the Gulf War began [5], In 1994 the 
United States Air Force admitted that some radars, including some mobile units, could 
detect the B-2 bomber [6], And as one would expect, the Russians are also developing 
their own anti-stealth technology [6], The bottom line is that if stealth platforms are to 
remain stealthy, they need to continue to decrease their observability. 
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The idea of stealth is comprised of several aspects of observability including 
infrared radiation, optical, acoustic, and radar echo. The objective is to minimize the 
observables associated with your platform so as to reduce the chances of being detected 
by the enemy. While current stealth platforms have incorporated reductions in each of 
these areas, the most emphasis has been placed on reduction of radar signature. This is a 
consequence of radar being the primary means of detection for the U.S. military, its allies, 
and its adversaries. That being the case, this thesis attempts to take one small step in the 
direction of radar cross section reduction by determining the induced source distribution 
on a body which produces scattered electromagnetic radiation. 

2. Goals of the Research 

Radar cross-section (RCS) reduction is a key element of low-observability. 

Current techniques used in determining the RCS of a platform rely on analysis of 
measurements made in the far-field. Generally these measurements provide a gross picture 
of the platform’s overall RCS as a function of viewing angle. This information enables 
engineers to then modify the design in an attempt to reduce its RCS. This becomes an 
iterative process with design modifications leading to measurements which lead to more 
design modifications. 

At the other end of the spectrum, measurements made using a probe on or near the 
outer surface of the body may influence the quantities being measured. In other words the 
act of measuring would induce errors into the quantities being sought. 

The focus of this study is to more accurately determine the source of scattering on 
a body. This will be done by evaluating the viability of back-propagating measurements 
to an axisymmetric body to image the source of scattering on the body. This analysis will 
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enable more precise location of scattering sources than does gross analysis of the far-field. 
Once the local scatterers have been identified it is then theoretically possible to remove or 
reduce them in the overall effort of radar cross-section reduction. 

The objective of this study is to investigate the viability of back-propagating 
electromagnetic field measurements to an axisymmetric body in order to determine the 
source distribution for radiated power from the body. If this is indeed possible, then the 
follow-on objectives include determining the range and spatial resolution at which 
measurements must be made in order to provide meaningful results. 



B. SOURCE IMAGING 

1. Acoustics Work with Cylinders 

A brief history of this subject starts with researchers attempting to determine 
sources of acoustic noise on a body. It was shown that acoustic intensity could be used to 
localize sources of sound on structures which radiate to the far field. A new quantity 
called the supersonic acoustic intensity vector was defined and its application to 
measurements on plate and cylinder-like structures was demonstrated [7], 

“Supersonic intensity is composed only of wave components which 
radiate to the far field (supersonic), with the non-radiating (subsonic) 
components eliminated. The normal component of this supersonic intensity 
vector, measured in the extreme near field or on the surface of the 
structure, provides an accurate tool for locating regions (“hot spots”) on 
the structure which radiate to the far field. Furthermore, the supersonic 
intensity provides an accurate quantification of these source regions, 
providing a ranking of the strength of the identified source regions as a 
function of frequency. This identification and ranking provides a powerful 
new tool in the understanding and control of radiated noise.” [7] 
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Coupled with the finding that acoustic source regions could be determined through 
back-propagation was that the results were obtained in two dimensions. Using a cylinder 
as the scattering body, scattering data measurements were made on an arc of a circle with 
diameter on the axis of the cylinder [8], In the spherical ( r,6,<j > ) coordinate system this 
translates to measurements over 0 at a constant <f > . The restricted case of axisymmetric 
fields and currents is the next level of difficulty before working with full three dimensional 
shapes. 

2. Extension to Electromagnetics for Cylinder Geometry Case 

The discoveries employing supersonic modes to enhance acoustic imaging on 
cylindrical shells has been translated into the realm of electromagnetics for the case of 
cylindrical geometry. It was shown that superluminal modes which satisfy |k 2 | < k and 
have faster-than-light axial propagation velocity provide time-average power to the far 
field. An optimal propagation deconvolution filter was implemented to remove the effects 
of the subluminal modes during back propagation. [9] 

3. Extension to General Non-separable Geometries Using Finite Element 

Methods 

The goal of this research is to extend the concepts developed in previously 
conducted research to bodies with non-separable geometries. By definition these bodies 
have shapes that do not lend themselves to explicit rigorous solutions using separable 
modes, such as cylindrical or spherical wave functions. A finite element method (FEM) is 
used to relate the conditions at the boundary (measuring region) to conditions on the 
surface of the body, in this case relating the scattered field to the source current that 
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generated it. The number of elements used is chosen so that the contour of the body is 
closely matched. Figure 1 shows a meridian plane of a sphere surrounded by a finite 
element mesh which connects the surface of the sphere to the arc of a circle with diameter 
on the axis of the sphere, as previously mentioned. In Figure 2, the body represented is a 
cone. The more contours on the body, the greater the number of elements required to 
closely match those contours. Ultimately the accuracy of the solution is determined by 
the number of elements used. 
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Figure 1 . Sphere (Meridian Plane) Surrounded by Finite Element Mesh 
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Figure 2. Cone (Meridian Plane) Surrounded by Finite Element Mesh 



In the chapters that follow, the source current on the surface of an axisymmetric 
body is determined through back-propagation of the scattered field it generates. In 
Chapter II, the scattered field (measured field) is determined through numerical integration 
of a surface current induced on an axisymmetric body by a locally placed elemental dipole. 
Chapter III details the back-propagation of the scattered field. Chapter IV provides an 
analysis of both the numerical integration results and the back-propagation results. 
Conclusions are presented in Chapter V. An appendix which includes all pertinent 
computer code is also included to assist in the understanding and continuation of this 
research. 
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n. GENERATION OF SCATTERED FIELD (MEASURED FIELD) 



A. BACKGROUND 

Before the scattered field can be back-propagated it must be measured. To 
simulate the process and control the sources of error, the measured field is computed to a 
specified level of accuracy. While reducing potential error associated with field 
measurements, generating fields leads to an additional hurdle, namely how to generate a 
scattered field from an arbitrary axisymmetric body. This problem is solved by inducing a 
surface current on the body and integrating it to find the resulting scattered field. 

Potential error is introduced as a result of integration accuracy. A sphere is first chosen as 
the arbitrary axisymmetric body in order to test the accuracy of the general source 
integration algorithm. 

A particular scattered field is generated by placing a radial directed elemental 
dipole in the vicinity of the metal sphere. The dipole induces a surface current J s on the 
sphere which in turn generates scattered fields E s and H s . The measured field H at 
some specified distance from the sphere is composed of the field scattered (H s ) from the 
sphere, and the field incident (H i ) from the elemental dipole. The ‘exact’ induced 
currents and scattered fields can be computed for this test case and used to validate the 

accuracy of the numerical integration and finite element algorithms. Calculation of the H 
field is considered in Section C. It is generated using the program DEPSPHR2.M found in 
the Appendix. Field integrations for axisymmetric currents are detailed in Section D and 
integration results are provided in Section E of this chapter. 
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B. SPECIALIZATION TO AXISYMMETRIC FIELD CASES 

Determining the source distribution on bodies with non-separable geometries is the 
ultimate objective of research in this area. This thesis reduces the problem from three 
dimensions to two dimensions by assuming axisymmetric fields generated by axisymmetric 
currents on a body of revolution. There are two special cases to be considered: the TE ^ 

case in which the E field and the surface current are transverse to $ , and its dual, the 
7M, case in which the H field is transverse to <f> and the surface current is <f> -directed. 
These cases are depicted in Figure 3 . 




Restricted Case: TE ^ 



Dual: TM^ 



Figure 3. Axisymmetric Field Cases 



This thesis will investigate the TE ^ case. The fields for this special case are. 



H(p,z)=H,(p,z)j 



(la) 



E(p,z) = E p (p,z)p + E z (p,z)z 



(lb) 
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The fields for the dual special case, TM + , are, 

E(p,z) = E 4 (p,z)<f> 



(2a) 



H(p,z) = H p (p,z)p + H 2 (p,z)z 



(2b) 



Maxwell’s curl equations are: (with a> = 2 zrf ) 



V x H = jcosE + J 



(3a) 



V x E = - jcopH 



(3b) 



Using Equation (3a) for the assumed field components in Equation (1) gives 




(4a) 



1 1 d 1 

j ( °E 2 = (p//,) —J z 

s p op r e 



(4b) 



These are the generating equations for E in terms of and known source current J t . 
Using Equation (3b) for the assumed fields in Equation (2) gives 



These are the generating equations for H in terms of E^ . 

C. OFFSET DIPOLE TO GENERATE SURFACE CURRENT ON SPHERE 
A radial directed dipole located in the vicinity of a metallic sphere is used to 
generate surface currents on the sphere. These surface currents are then integrated to find 
the scattered field at a specified distance. The dipole case is also used to test the 




(5a) 




(5b) 
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integration accuracy by providing a nearly exact scattered field. Figure 4 shows an 
elemental dipole located in the vicinity of a metal sphere in the spherical coordinate 
system. 




Figure 4. Dipole in Vicinity of Metal Sphere 

The procedure for determining the surface current consists of several steps. The 
E t and H i fields due to the dipole are represented in ( r,0,<f > ) coordinates of the centered 
sphere. The scattered fields E s and H s due to the induced currents on the sphere are 
then expanded using spherical harmonic expansions with unknown coefficients. Enforcing 



10 



fx(£,. + E S ) = 0 => Zsun = 0 on the sphere surface allows solution for the expansion 
coefficients. Finally J s =rx (//. +H s ) is used to find induced surface currents which 

produce E s and H s fields scattered from the sphere. 

An elemental dipole placed at the origin produces the following fields [10]: 



zr a 

E " =~ cos *- 



_ 7.W/ . . 

E >. = T^~ suie ‘ 



1 



J 



fir, 



, - iP'd 



JJL + J_ 



fir c 



jf> r d 



Idl 

H 4 = sin 

** 4 n 



J£ + ± 

r} 



jfiU 



(6a) 



(6b) 



(6c) 



The position of the dipole must then be translated to the position z = z 0 . This 
translation is accomplished using the following transformation formulas: 



r d = yl r * +z o ~ 2rz o cos & 


(7a) 


r cos 6 - z n 

COS0, - 

r a 


(7b) 


. rsin<? 

sin#, - 

r* 


(7c) 


Eg = Eg d cos (0 -0 d )~ E rj sin(6» - Q d ) 


(8a) 


E r = Eg d s\n(6-0 d ) + E, d cos((9 - 6 d ) 


(8b) 


ii 


(8c) 
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The transformation formulas in Equations (7) and (8) are then used with 
Equation (6) to express the fields due to the dipole on the surface of the metal sphere, 



E gj (a,0 ) , E ri (a,0 ) and H^{a,6 ) . This field will induce a surface current J = J g {6)6 
on the surface of the sphere which will, in turn, generate scattered fields E gs , E rs , and 
Hj s . The tangential E field at the sphere surface will exactly cancel the tangential E 
field due to the dipole thus producing a total tangential E field of zero 



This J g {6) is what needs to be integrated to find the radiated H^ s (r,0) at locations off 
the surface of the metal sphere. The result is H^ s which will be used to compare results of 

the integration covered in the following section. 

Note that although J e is produced by the total on the metal surface, it 

generates only the scattered field. The scattered field has an E 6s which exactly cancels 
out the &0i of the incident dipole field on the sphere surface. 

To determine the “exact” solution, the scattered fields can be expressed as 
weighted sums of spherical harmonics [11], 



E e ,(a,V) = -E„S a ,6) 



(9a) 



The surface current will also satisfy J = h x H, ota i on the surface: 





(9b) 




(10a) 
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H t ,(r,9) = costf) 



(10b) 



where H^ 2> (J3r) = J„(J3r)~ jY n (J3r) is the “Riccati” spherical Hankel function of the 
second kind representing outbound waves. P n ‘(cos^) is the “associated Legendre 
function” and is the m = 1 case of P n m (cos^) . To numerically solve the problem, 
truncate the series in Equation (10) and solve for N values of complex a n . Using 

Equation (9a) and defining D n = gives 

(3a 



N 



j J hH a n D n P n(cos0) = -E gi (a,0) 



( 11 ) 



n=\ 



The a n are found by use of orthogonality of the associated Legendre functions, 
P n m {cos0 ) . For the special case m = 1 , 



A 

j P( (cos#)/} 1 (cos#) sin Odd = 



2n(n + 1) 



2w + l 
0 l*n 



l = n 



( 12 ) 



We can sift out coefficients: 



JnP,D m = - f E gi ( a,0)P' n (cos #) sin 6d9 

2/7 + 1 i 



(13a) 



a = 



JHfi, 2n(n + 1) 



(13b) 



The I n integration can be performed numerically to any accuracy since E di (a,6) and 



P n '(cos#) can be obtained using as many 0 -points as desired. Once the a n are found the 
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J e {0) can be computed at any 0 value by use of Equation (9b) with H+ s given by 
Equation (10b). 

D. FIELD INTEGRATION FOR AXISYMMETRIC CURRENTS 

As noted in the introduction, the first step is to measure the field that will be back- 
propagated to the surface of the axisymmetric body. For our purposes we will assume an 
axisymmetric surface current distribution on the body. The “measured” fields will actually 
be determined through the integration of the surface current. In this context the scattered 
field is generated through integration rather than through illumination and measurement. 

The surface of an arbitrary axisymmetric body of revolution can be defined by use 
of the generating contour p'(z') using the ( p,<f>,z ) cylindrical coordinate system as 
shown in Figure 5. Source points on the body are defined by individual (p,z) pairs. The 
field point in space is located at (fi,<j>,z) . The axisymmetric surface current distribution is 

A A A A 

J s (r') = JtiP't 2 ')* + J • The unit vectors t and <f> are tangent to the surface at 
each point and both J , and J + are constant with changing <j > . 
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t 

Figure 5. Axisymmetric Body of Revolution 



The fields will be numerically evaluated via the vector potential formulation given 



by [11]: 



H(r) = — V x A(r) 
Mo 



(14a) 



E - — - — V xH = jo)A - — V(V • A) 
jcos 0 co^ie 



(14b) 



where 



JJ J ‘ ( - r "> 



— — e 



-m 



surface 



R 



-dS' 



(15) 
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with R-r-r' and R = |r - r '\ . For field points and source points with coordinates 
(fi,<j>,z) and (j)' ,(f>' ,z') the law of cosines gives 

R = ^/p 2 + p' 2 - 2 pp' cos(^ - ^') + (z - z') 2 (16) 

The curl operator in Equation (14a) is taken inside the integral in Equation (15) and note 
that it only operates on the e ^/r term (unprimed position in R = |r - r'| ) 

= jj vf— -^-J x~J s (p)dS' (17a) 

surface ^ ^ 

where 



V 




= -R 



1 + 7 /^ 



R 



ipR 



(17b) 



Thus the general solution is: 



»('>= 47 



surface 




\ 1 1 ££ \ 

R 3 ) 



e~ ipR dS ' 



(18) 



The next step is to derive the explicit integrations to be performed. To do so we 
select r in the x-z plane, with <f> = 0 as shown in Figure 6. Since J s and the resulting 
fields so derived are axisymmetric (not functions of ^ ) we lose no generality by doing this. 
In fact, if we find H x ,H y and H z at this field point we note that for other locations where 
<f> -*■ 0 the substitutions H x —> H p ;H y —> and H 2 —> H : apply. Also note that 



r - px + zz (19a) 

r' = p' cos<f>'x + p' s\n<f>'y + z'z (19b) 
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R - r - r' - (p-p'cos^')x-p' sin^'_y + (z -z')z 



(19c) 



R = |i?| = -^/p 2 + p' 2 - 2 ~pp' cos^' + (r - z') 2 (19d) 




Figure 6. Geometry for Explicit Integrations 



The field integrations will be performed numerically by breaking the axisymmetric 
surface into a large number of flat circular rings as shown in Figure 7a. These rings are 
stacked to approximate the surface of revolution. Each ring is “flat” in the direction of 

t and curved in the direction of $ as depicted in Figure 7b. 
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Figure 7a: Stacked Circular Rings form Axisymmetric Body 



z 




Figure 7. Decomposition of Axisymmetric Body 
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Consider now evaluating H(r ) in Equation (18) due only to J s (r') = J t (s)t over 
a single ring as in Figure 7b. Over the ring we will approximate J t (s ) as piecewise linear, 



specified using its values at the top and bottom of the ring, J t ( 5 , ) and J t (s 2 ) via 

Is - s . ) (s - 5 . ) 

\ S l S 2 ) \ 2 ^ 1 ) 



( 20 ) 



This equation is a straight line interpolating function between known values at s , ands 2 . 

A single ring is shown in Figure 8a where s is the path length over the ringed surface from 
the + z axis. The piecewise linear variation of J, (s ) between s ] and s 2 is shown in 

Figure 8b. 

To evaluate Equation (18) for J s ( r> ) = *4 (s)t over a single ring, we note that 
t = cosacos^'x + cosasin^'_y-sina£ (21) 

where cos a = tp and - sin a = t • z . After simplifying, 

txR = x(t y R, - t z R y ) + y(t 2 R x - t x R,) + z{t x R y - t y R x ) (22) 

and substituting Equation (22) into Equation (18) with restriction to one ring gives, 



1 * 2 

H x = — | J t (s)[(z - z’) cos a - p’ sin a]F s (s)p'ds 



(23a) 



1 ^ 

H y = — J J t (^){[p' sin a - (z - z ') cosa]F c (s') - p sin aF } (s)^p'ds (23b) 



H z =-^-\j t (s)p cos aF s (s)p’ds 



An 



(23c) 
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Figure 8b: Piecewise Linear Variation J t (s) 
Figure 8. Surface Current Variation Over One Ring 



where 



* 

F s( s )= J sin ^ 



R 3 



F c (s) = ] cos<p'[-- + ^ R ^ e 



j e~ J/ *dp 


(24a) 


\- jf<R d<f)' 


(24b) 
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(24c) 




As we vary <)>' while holding 5 constant (e.g. p',z' constant), R{<j>') is an even 
function: R(-<j > ') = R(<fi') • Therefore ^ ^ is also an even function of <}>' . 



Since sin^' is an odd function of <j>' the integrand in Equation (24a) is also odd giving 
F s (s) = 0. We are left with only H in (10b) since H x = H 2 = 0 . 



Extending this to any position r , we have for J s = J t {s)t , 






(25a) 



where 



=-^— Jj J (5){[/j'sina-(z-z')cosa]/ r c (5)-psinaF I (5 , )J/7'<iy (25b) 



and 



Fc(s) 



= 2|cosf(-^ 



m 



(26a) 






(26b) 



E. CENTERED SPHERE INTEGRATION TEST CASES 

A centered sphere was chosen as the test case to determine the accuracy of the 
integration algorithm. The integration algorithm was developed using the trapezoidal rule 
Nine cases were tested with sphere sizes ranging from a radius of one wavelength to ten 

wavelengths. For each sphere size the H field was generated at three different locations: 
0.2, 1 and 3 wavelengths from the outer surface of the sphere. 
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The offset distance of the dipole was kept constant for each of the three cases associated 
with the three different sphere sizes. The number of modes used was determined by the 
following: 

N=integer(2 ka + 2), (27) 

where k is the wave number (2 n / X) and a is the sphere radius. In each of the nine 
cases the frequency of operation was 300MHz. The RMS error was determined by 

comparing the H field determined through integration to the H field exact solution found 
using the procedure developed in Section II.C. Each of the cases was tested with 
increasingly fine segmentation in 0 and <f> on the sphere surface until the RMS error was 
in the neighborhood of one percent. In order to achieve the desired integration accuracy, 
the number of integration points on the body was increased. This relates to the spatial 

resolution required when measuring the H field in order to achieve acceptable results 
when back propagating that field. As the number of integration points on the body 

increases, the measured H field becomes more accurate. The more accurate the 
measured field, the better the surface current on the sphere can be resolved. The test 
cases summarized in Table 1 are chosen based solely on the RMS error associated with 
each being in the neighborhood of one percent. 

The results shown in Table 1 indicate the RMS error for source integrated//^, 

observed at the “field distance” radius. The RMS error is a function of the number of 
surface points at which the surface current was integrated. Figures 9 through 12 show 
how the integration results converge for several of the cases shown in Table 1. For the 
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CENTERED SPHERE INTEGRATION RESULTS 




Case 


Sphere 

Radius 


Field 

Distance 

(*) 


<f> 

Segments 


e 

Segments 


Dipole 

Offset 

W 


Modes 


RMS Error 
(%) 


la 


1 


1.2 


21 


21 


5 


14 


1.1490 


lb 


1 


2 


21 


21 


5 


14 


1.1620 


lc 


1 


4 


21 


21 


5 


14 


1.1510 


2a 


5 


5.2 


78 


78 


9 


64 


1.0310 


2b 


5 


6 


78 


78 


9 


64 


0.1220 


2c 


5 


8 


78 


78 


9 


64 


0.1250 


3a 


10 


10.2 


157 


157 


14 


128 


0.8268 


3b 


10 


11 


157 


157 


14 


128 


0.0456 


3c 


10 


13 


157 


157 


14 


128 


0.0459 



Table 1. Centered Sphere Integration Results 



cases in Table 1, the number of <f) and 6 segments is based on the radius of the sphere. 
However, the number of <j> segments remains constant as the radius of each ring (see 
Figure 7) decreases near both poles of the sphere. This finer (j> segmentation gives 
increasingly narrower patches of surface area over which the integration is performed. 
Consequently, the RMS error does not appear to decrease significantly as the number of 
<j> segments is increased. This is a result of the increased <j> segmentation that is “built in” 
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as the radius of the rings decrease. An increase in RMS error is observed for the 
decreased field distance from the sphere in Cases 2a and 3 a. Since the field distance is 
only 0.2/1 from the sphere surface in these cases, the large relative contribution and rapid 
variation of inverse distance terms in the integrand of Equation (18) for surface points 
closest to the field point require increasingly fine segmentation for accurate numerical 
integration. 



Integration Convergence: Casela 
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Hphi RMS Error Hphj RMS Error 



Integration Convergence: Caselb 




10 



Integration Convergence: Case2a 




70 
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Integration Convergence: Case2b 




The magnitude and phase of the H field at the measured distance are presented in 

Figures 13 through 30. These figures help to tell the complete story of the computed H 
field. In each figure the “exact” quantity computed via series solution is depicted by a line 
and the “measured” quantity is depicted by dots. Note that in both the magnitude and the 
phase plots the measured quantity closely matches the exact quantity. These results 

indicate that the integration produces an H field with minimal error. This knowledge will 
help isolate the source of error that may present itself due to the back-propagation 
process, which will be described in Chapter III. 
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Hexact (line) & Hinteg (dots): Casela; RMS Error= 1.149% 




Hexact (line) & Hinteg (dots): Casela 
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Theta (degrees) 



Figure 15. Magnitude Comparison: Case lb 



Hexact (line) & Hinteg (dots): Caselb 
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Hexact (line) & Hinteg (dots): Caselc; RMS Error= 1.151% 




Hexact (line) & Hinteg (dots): Caselc 
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Hexact (line) & Hinteg (dots): Case2a; RMS Error= 1.031% 




Hexact (line) & Hinteg (dots): Case2a 




30 




Phase Angle (radians) 



Hexact (line) & Hinteg (dots): Case2b; RMS Error= 0.122% 




Hexact (line) & Hinteg (dots): Case2b 
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Hexact (line) & Hinteg (dots): Case2c; RMS Error= 0.1259% 




Hexact (line) & Hinteg (dots): Case2c 
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Hexact (line) & Hinteg (dots): Case3a; RMS Error= 0.8268% 




Hexact (line) & Hinteg (dots): Case3a 
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Hexact (line) & Hinteg (dots): Case3b; RMS Error= 0.04562% 




Hexact (line) & Hinteg (dots): Case3b 
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Hexact (line) & Hinteg (dots): Case3c; RMS Error= 0.04596% 




Hexact (line) & Hinteg (dots): Case3c 
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m. BACK PROPAGATION OF SCATTERED FIELD (MEASURED FIELD) 



A. BACKGROUND 

A finite element method is used to back-propagate the measured field to the 
surface of the sphere. The measured fields that are back-propagated are those determined 
in Chapter II and summarized in Table 1. Once the measured field is back-propagated, the 
surface current is determined using Equation (9b). The surface current due to back- 
propagation (Jg^ ) is compared to the original surface current (Jo ^, ) due to the 
elemental dipole which is computed using a spherical harmonic series, per Section II. C. 
The following section describes the finite element formulation. 

B. FINITE ELEMENT FORMULATION FOR AXI SYMMETRIC FIELDS 

For the special case being considered, with no <j) variation in materials or sources, 
and J j - 0 , the fields are derived in Section II. B., results of which are repeated here for 
convienence. 



H{p,z) = H,(p,z)<f> 



(28a) 



E{p,z) = E p (p,z)p + E 2 (p,z)z 



(28b) 



Maxwell’s curl equations are: (with co = 2jrf ) 



V x H = jcosE + J 



(29a) 



V x E = -jatfiH 



(29b) 
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Using Equation (29a) for assumed fields in Equation (28) gives 



• p 1 ; 

; “ E ' = 7T'7' / ' 

1 1 1 
]aE ^l-pTp ipH > ) -~s J ’ 



(30a) 

(30b) 



These are the generating equations for E in terms of and known source 
current J . In this case Equation (29b) simplifies to: 



jo)(V xE)-</> = ja> 



cE p cE z 



dz dp 



= co 2 pH+ 



(31) 



Substituting the E fields from Equation (30) into Equation (31) gives the partial 
differential equation (PDE) satisfied by : 



d 


" 1 d , T ^ 


d 


\\cHA 


dp 


IwH 


^ dz 


s dz 



d ,1 



d .1 






(32) 



dp 8 2 dz e 
where e - s(p,z ) is allowed to be to be inhomogeneous. 

The PDE in Equation (32) will be numerically solved using the finite element 
method (FEM), where H t is approximated with electrically small triangular element 

regions in the ( p,z ) plane. [12] 

A variational Euler-Lagrange algorithm is used to implement the finite element 
method. This algorithm replaces the solution of the PDE in Equation (32) with finding the 
stationary point in complex function space of a quadratic (in H , ) functional. We first 

expand the field using pyramidal basis functions, u n (p,z), defined for each node in the 
mesh and having support of the surrounding triangular elements. 



38 



( 33 ) 



H+{p,z) = &A(P-z) 

n= 1 

and then substitute Equation (33) into the functional. The stationary solution is found by 
differentiating with respect to the unknown values of H ^ , located at the m-th nodes, to 

yield the linear system 

N re 1 

2X j] - V(pw m ) • VO„ ) - /: 2 pu m u n dpdz = 0 (34) 

n ~ 1 overlapping elements • 

Defining the double integrals indexed on m,n as T (m,n) , the linear system in Equation 
(34) can be rewritten as, 

Z hj(m,n) = - ^h n T(m,n) (35) 

n inside boundary n on boundary 



->na(n) hb(n ) 



0“ 


~K ' 






A 


h 


= 


B 


_0 


1 

| 
J 




- 



— J U- -rJ -J 

Square N m xN s 

( ve/y sparse ) ( sparse ) 



boundary 

values 



(36) 



The m-th row of A and B is formed by Equation (35) with the n-th term 
corresponding to a column of A if h n is unknown. The appropriate column of B is filled 
if h n is a boundary node. There is no contribution to A or B for nodes on the z-axis 
{p = 0) since H^(p = 0) = 0 and those nodal values of h n are known. [13] 



C. FINITE ELEMENT METHOD LIMITATIONS 

The finite element method was successful for back-propagation of fields measured 
at distances less than 0.25 wavelengths from the surface of the sphere. When the field was 
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measured at distances greater than 0.25 wavelengths, errors due to back-propagation 
increased significantly. The same cases for which the measured field was generated were 
used for back-propagation. These results are summarized in Table 2. The program 
FEM2.M was used to determine the surface current due to back-propagation. It can be 
found in the Appendix. 



BACK-PROPAGATION RESULTS. FINITE ELEMENT METHOD 


Case 


Sphere 


Field 


Modes 


H* 


J* 




Radius 


Distance 




RMS Error 


RMS Error 




<*) 


(*) 




(%) 


(%) 


FEMla 


1 


1.2 


14 


0.0590 


0.9976 


FEMlb 


1 


2 


24 


0.0338 


30.09 


FEMlc 


1 


4 


50 


0.0147 


25.87 


FEM2a 


5 


5.2 


64 


0.0043 


0.7406 


FEM2b 


5 


6 


74 


0.0052 


57.18 


FEM2d 


5 


5.26 


64 


0.0042 


662.9 


FEM3a 


10 


10.2 


128 


0.0017 


0.7333 


FEM3d 


10 


10.26 


128 


0.0063 


5.547 



Table 2. Centered Sphere Back-Propagation Results: Finite Element Method 
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The magnitude and phase of the surface current on the PEC sphere, as determined 
through the finite element method, is shown in Figures 3 1 through 42. These figures help 
to tell the complete story of the accuracy in back propagating to obtain the current. In 
each figure the “exact” quantity is depicted by a solid line and the calculated quantity is 
depicted by dots. Note that in both the magnitude and phase plots the calculated quantity 
closely matches the exact quantity for the cases in which the H field is measured less than 
0.25 wavelengths away from the surface of the sphere. In the cases where the distance is 
one wavelength, the error is immense. Cases FEM2d and FEM3d have been added to the 
original nine test cases in order to demonstrate what happens when the measurement 
distance is increased to 0.26 wavelengths from the surface of the sphere. The following 
section will investigate the cause of this error. 
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Jexact (line) & Jback (dots): FEMIa; RMS Error= 0.9976% 




Jexact (line) & Jback (dots): FEMIa 
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Phase Angle (radians) 



Jexact (line) & Jback (dots): FEMIb; RMS Error= 30.09% 




Jexact (line) & Jback (dots): FEMIb 




Figure 34. Phase Comparison: FEMIb 
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Phase Angle (radians) 



Jexact (line) & Jback (dots): FEMIc; RMS Error= 25.87% 




Jexact (line) & Jback (dots): FEMIc 




Figure 36. Phase Comparison: FEMIc 
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Jexact (line) & Jback (dots): FEM2a; RMS Error= 0.7406% 




Jexact (line) & Jback (dots): FEM2a 




Figure 38. Phase Comparison: FEM2a 
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Jexact (line) & Jback (dots): FEM2b; RMS Error= 57.18% 




Jexact (line) & Jback (dots): FEM2b 
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Jexact (line) & Jback (dots): FEM2d; RMS Error= 662.9% 




Jexact (line) & Jback (dots): FEM2d 
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Jexact (line) & Jback (dots): FEM3a; RMS Error= 0.7333% 




Jexact (line) & Jback (dots): FEM3a 
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Jexact (line) & Jback (dots): FEM3d; RMS Error= 5.547% 




Figure 45. Magnitude Comparison: FEM3d 



Jexact (line) & Jback (dots): FEM3d 
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When back-propagating measurements are made at distances greater than 0.25 
wavelengths from the surface of the sphere using the FEM approach, the RMS error of the 
predicted surface current tends to increase dramatically relative to the situation of smaller 
distances. This appears to result from resonant modes injecting themselves into the 
solution in varying degrees. Figure 47 shows the single mode RMS error for a centered 
PEC sphere used in case FEM la. Note that the RMS error associated with each mode is 
relatively low; all but two are below one percent. The results shown Table 2 indicate that 
back-propagation gave a solution with less than one percent RMS error. Figure 48 shows 
the single mode RMS error for case FEMlb. In this case the measured field was one 
wavelength from the surface of the sphere. Note that the lower order modes generally 
have RMS error in the neighborhood of five percent or less, except for mode n = 5 which 
is nearly fifty percent RMS error. This error “spike” indicates a resonant mode that 
significantly alters the solution. As can be seen in Figures 48 through 53, the resonant 
modes seem to appear in only the cases where the measurement distance is greater than 
0.25 wavelengths from the surface of the sphere. 
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Percent 



Single Mode RMS Error for Centered PEC Sphere: FEMIa 




Mode Index "n" 



Figure 47. Single Mode RMS Error: Case FEMIa 



Single Mode RMS Error for Centered PEC Sphere: FEMIb 




Figure 48. Single Mode RMS Error: Case FEMIb 
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Percent 



Single Mode RMS Error for Centered PEC Sphere: FEM2a 




Figure 49. Single Mode RMS Error: Case FEM2a 



Single Mode RMS Error for Centered PEC Sphere: FEM2b 
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Percent 



Single Mode RMS Error for Centered PEC Sphere: FEM2d 




Single Mode RMS Error for Centered PEC Sphere: FEM3a 




Figure 52. Single Mode RMS Error: Case FEM3a 
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Single Mode RMS Error for Centered PEC Sphere: FEM3d 




Another way to view this phenomenon is by looking at the product solution parts 
of H^{r,G) = g(r) ■ h(0) for a given mode, n. The radial variation, g(r ) , to be shown in 

the upper subplots, indicates the contribution at the PEC sphere surface due to a constant 
contribution of g(b) = 1 on the outer boundary. Figures 54 and 55 show the radial 
variation for modes n = 4 and n = 5, respectively, for the case FEMlb. Note that mode 4 
provides relatively small contribution at the PEC sphere surface. Mode 5, however, 
indicates an extremely large contribution at the PEC sphere surface. Two additional cases 
have been added here to further illustrate this observation. Cases FEM2d and FEM3d 
shown in Figures 56 and 57, respectively, are both cases in which the measured field is 
0.26 wavelengths from the surface of the sphere. Note the large contribution at the 
surface of the sphere due to the resonant modes captured in each of these figures. 
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Theta Variation , Radial Var,ation Theta Variation Radial Variation 



Product Solution Parts of Hphi(r.theta). FEMIb; n=4 




Theta (degrees) 

Figure 54. Product Solution Parts: Case FEMIb, n = 4 



Product Solution Parts of Hphi(r, theta): FEMIb; n=5 




[I . ' 

0 50 100 150 



Theta (degrees) 

Figure 55. Product Solution Parts: Case FEMIb, n = 5 
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Product Solution Parts of Hphi(r, theta): FEM2d; n=8 




Figure 56. Product Solution Parts: Case FEM2b, n = 8 



Product Solution Parts of Hphi(r, theta): FEM3d; n=16 





Figure 57. Product Solution Parts: Case FEM3b, n = 16 
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Additional insight is provided by viewing the resonant modes that appear when the 
field is measured at varying distances from a fixed size sphere. Figure 58 considers modes 
n = 1 to 100. Each line in Figure 58 displays the contribution of g(a ) from a single mode 
at the sphere surface as a result of a constant contribution of g(b) = 1 at the boundary 
distance. The “best radius,” b = 5.641 \ in this case, provides the minimum g(a) from 
all modes applied as g(b) = 1 . From the perspective of minimal resonance effects, this is 
the best location to measure the field. 




Figure 58. Resonant Effects at Various Outer Boundary Locations 

D. TRANSFER MATRIX FOR CENTERED SPHERE 

As can be seen in Table 2 the results due to back-propagation are poor for 
measurements made greater than approximately 0.25 wavelengths from the surface of the 
sphere. A transfer matrix using spherical harmonics for a centered sphere was developed 
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to accurately back propagate the field on the outer boundary to that on the sphere, 
without the bothersome injection of resonant modes as appear with the FEM solution. 

This method will give some indication as to the spatial resolution required when measuring 
the field so that the source current on the sphere can be accurately determined. 

Figure 58 shows a meridian plane for a centered sphere of radius a enclosed by an 

outer boundary arc of radius b . For any radius r , where a < r < b , the total H field can 
be determined by: [14] 



N 






J n (kr) H™(kr)~ , 
nV - + b_ V P„ (cos 0) 






kr 



kr 



(37) 



where: 



J'Skd) 



forces Eg - 0 for each mode on the surface of the sphere. 



(38) 



A specified / ( 0 ) = H^(b,0) will provide a set . The a n coefficients are 

determined by moment matching, using orthogonality of Legendre functions. This 

technique is detailed in [14 ]. The transfer matrix H s is developed such that 

H+(a,Q q ) = H s (g,p)- Hf(b,Q p ) (39) 
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Centered Sphere radius = 2; Arc radius = 5 




Figure 59. Meridian Plane of Centered Sphere 

where H^{b,9 p ) is a column array of the 6 p sampled field at radius b , as determined in 
Chapter II. The array product results in the column array of H* at radius a which 

corresponds to the surface of the sphere. From this result the surface current is then 
determined using Equation (9b). The test cases using the transfer matrix were generated 
using the program HSPHERE1.M which can be found in the Appendix. Table 3 tabulates 
the results and compares them to those obtained using the finite element method. It 
should be noted that the transfer matrix works only for a centered sphere and can not be 
applied to any arbitrary axisymmetric body. 
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SURFACE CURRENT RESULTS USING TRANSFER MATRIX 


Case 


Sphere 


Field 


Modes 


6 


H* 


U 




Radius 


Distance 




Segments 


RMS Error 


RMS Error 




(*) 


<*) 






(%) 


(%) 


la 


1 


1.2 


14 


94 


0.0590 


0.8085 


lb 


1 


2 


24 


125 


0.0338 


0.6271 


1c 


1 


4 


50 


188 


0.0147 


0.7349 


2a 


5 


5.2 


64 


408 


0.0043 


0.7946 


2b 


5 


6 


74 


376 


0.0052 


0.9453 


2d 


5 


5.26 


64 


408 


0.0042 


0.7764 


3a 


10 


10.2 


128 


801 


0.0017 


0.7885 


3d 


10 


10.26 


128 


801 


0.0063 


0.7827 



Table 3. Surface Current Results Using Transfer Matrix 



The magnitude and phase of the surface current on the PEC sphere, as determined 
by the transfer array, is shown in Figures 60 through 77. Once again these figures help to 
tell the complete story of the surface current. In each figure the exact quantity is depicted 
by a solid line and the calculated quantity is depicted by dots. Note that in both the 
magnitude and phase plots the calculated quantity closely matches the exact quantity. 
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Jexact (line) & Jback (dots): Casela; RMS Error= 0.8085% 




Jexact (line) & Jback (dots): Casela 
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Jexact (line) & Jback (dots): Caselb; RMS Error= 0.6271% 




Jexact (line) & Jback (dots): Caselb 
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Jexact (line) & Jback (dots): Caselc; RMS Error= 0.7349% 




Jexact (line) & Jback (dots): Caselc 
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Jexact (line) & Jback (dots): Case2a; RMS Error= 0.7945% 




Jexact (line) & Jback (dots): Case2a 
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Jexact (line) & Jback (dots): Case2b; RMS Error= 0.9453% 




Jexact (line) & Jback (dots): Case2b 
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Jexact (line) & Jback (dots): Case2c; RMS Error= 0.5991% 




Jexact (line) & Jback (dots): Case2c 
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Jexact (line) & Jback (dots): Case3a; RMS Error= 0.7885% 




Jexact (line) & Jback (dots): Case3a 
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Theta (degrees) 



Figure 74. Magnitude Comparison: Case3b 



Jexact (line) & Jback (dots): Case3b 
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Jexact (line) & Jback (dots): Case3c; RMS Error= 0.7846% 




Figure 76. Magnitude Comparison: Case3c 



Jexact (line) & Jback (dots): Case3c 
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IV. ANALYSIS OF RESULTS 



A. THE INTEGRATION ALGORITHM 

Since the “measured” fields to be back-propagated were to be determined through 
integration, the first part of this thesis was dedicated to generating accurate fields through 
numerical integration. As is shown in Table 1, properly performed numerical integration 
provides highly accurate measured fields. Figures 9 through 12 provide several examples 
of how quickly the integration converges. In each case the percent RMS error is 
determined essentially by the number of points on the sphere surface over which the 
integration takes place. For the cases included in this thesis the number of points was 
chosen such that the error due to back-propagation, as determined by the transfer matrix, 
would be kept below one percent. 

Figures 13 through 30 show the magnitude and phase of the measured field 
determined through numerical integration. For each figure the measured field is over-laid 
on the exact field. Visual inspection indicates the magnitudes are nearly the same over the 
entire range of theta. RMS error computations indicate a difference of only a few 
hundredths of one percent. Visual inspection of the phase comparison gives similar results 
- both the computed and the exact are nearly the same. 

As the number of surface points increase for a constant sphere size, the accuracy 

of the measured H field increases. This result takes the form of decreasing RMS error as 
shown in the integration convergence of Figures 9 through 12. As the measurement 
distance increases for a constant sphere size, the number of points required for accurate 
integration results decreases. This can be seen by comparing Cases 2a and 2b shown in 
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Figures 1 1 and 12. For the larger measurement distance of Case 2b the RMS error is 
significantly smaller using the same number of points as Case2a. 

B. FINITE ELEMENT METHOD PROBLEMS 

The finite element method was used successfully to determine the source current 
on the PEC sphere when the field was measured at distances less than 0.25 wavelengths 
from the sphere surface. Table 2 summarizes the test cases used in this thesis. 

The large errors that arise when the field to be back-propagated was measured at 
distances greater than 0.25 wavelengths from the sphere surface appear to be caused by 
resonant modes. These modes are shown for several cases in Figures 48, 50, 52 and 53. 
When back-propagated, these modes provide overwhelmingly disproportionate 
contributions at the sphere surface. These contributions are shown for a few 
representative cases in Figures 55, 56, and 57. It appears these resonant modes come 
about as a result of the boundary conditions placed on the outer boundary, the location at 
which the field is measured. As can be seen in Figures 48, 50, 51 and 53, these resonant 
modes appear as sharp spikes in the plot. 

C. TRANSFER MATRIX 

The transfer matrix that was used to determine the source current on the body 
gave results with less than one percent RMS error. The low error is a result of the number 
of surface points used to determine the measured field, and the number of points at which 
the field was measured at the outer boundary. The relationship between number of 
surface points and numerical integration accuracy was discussed in Section A. As can be 
seen in Equation (18), the measured field varies inversely with distance from segment to 
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field point. This inverse distance variation is the main point of concern when R is small at 
the closest point of approach. There is a large relative contribution to the integration 
from this surface region and the variation in the field is rapid. The number of segments on 
the sphere surface in this near-region must be increased to ensure the accuracy of 
numerical integration. 
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V. 



CONCLUSIONS 



The objective of this research was to investigate the viability of back-propagating 
electromagnetic field measurements to an axisymmetric body in order to determine the 
source distribution for radiated power from the body. The method for achieving this 
objective was straightforward. First, the “measured fields” were simulated through 
numerical integration. Second, these fields were input into a finite element algorithm to 
determine the source current on a PEC axisymmetric body. The original thrust of the 
thesis was to investigate the ranges and resolution at which the fields could be measured in 
order to provide a usable level of accuracy in determining the source distribution on the 
body. The test cases were to include various “arbitrary” axisymmetric shapes to include 
cones, offset cones, offset spheres, and various cone-cylinder-sphere combinations. The 
centered sphere was to be used for very basic test purposes only. 

Along the way to the objective it was found that the finite element method 
provided inaccurate results in the centered sphere tests for all but very limited cases. This 
discovery prompted a redirection of efforts to investigate the sources of error. 
Consequently, progress towards the original goals was sidetracked; but, that is often the 
nature of true research into the unknown. 

The integration program developed to determine the measured fields was found to 
converge to an accurate solution. Although the shape tested was a sphere, the theory 
developed in Chapter II using stacked circular rings can be easily extended to arbitrary 
axisymmetric shapes. 

The finite element method was selected as a means to determine the source current 
on the sphere through back-propagation. The errors associated with this method for the 
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centered sphere were unexpected. Consequently, solving these associated resonance 
problems was not intended to be the focus of this thesis. Further testing of the FEM 
solution will be done as part of extensions to this effort in follow-on thesis efforts and 
techniques for mitigating resonance effects will be explored. 
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APPENDIX - MATLAB CODES 



%thesis2.M by Dan Wawrzyniak 5/10/97 
% updated 5/14/97 
% Combines all programs into one 
% Includes checking integration error first 
% Allows user to save data in a .mat file for future use 
clear all; 

case = input('What case are you running? 
fO = input('Enter frequency in MHz: '); 
lambda = 3e2/f0; 

disp(['Wavelength = ',num2str(lambda),' meters']); 

al = input('Enter metal sphere radius in terms of wavelength: ');a = al*lambda; 

rl = input(Enter outer boundary in terms of wavelength: ');r = r 1 *lambda; 

eta = 120*pi; % free-space Z_0 

Idl =1; % dipole moment 

B = 2* pi/lambda; % wavenumber 

zOl = input(Enter z-position of elemental dipole in terms of wavelength: '); 
zO =z01*lambda; 

Nthl = input('Enter number of theta points: '); Nth = Nth 1-1; 
ppwl = inputCEnter number points per wavelength: '); 
dth=pi/Nth; th=(0:dth:pi)'; thdeg=180*th/pi; 
sth=sin(th); cth=cos(th); 

% computing rd, sin(thd), and cos(thd) for sphere surface (r=a) 
rd=sqrt(a*a+z0*z0-2*a*z0*cth); sthd=a*sth./rd; cthd=(a*cth-zO)./rd; 

% computing field components in dipole centered coordinates 
Fth=Idl*exp(-j*B*rd)./(4*pi*rd); 

Hpi=Fth.*sthd.*(j*B + l./rd); 

Etd=eta*Fth.*sthd.*(j*B + l./rd - j./(B*rd.*rd)); 

Erd=2*eta*Fth.*cthd.*(l./rd - j./(B*rd.*rd)); 

% translating incident field spherical components 
cdth=cth.*cthd + sth.*sthd; sdth=sth.*cthd - cth.*sthd; 

Eti=Etd.*cdth - Erd.*sdth; % note that Hpi=Hpd 
% Computing spherical harmonic coefficients for scattered field 
% Ets expansion which minimizes the total sphere surface 
% tangential field = Eti+Ets 
Nmax = fix(2*B*r); 
disp([Nmax = ’,num2str(Nmax)]); 

N = input('Enter number of spherical harmonics: '); 

Sn=-dth*Eti.*sth; Pnl=legpol2(cth,N); In=(Pnl.')*Sn; 
n=(l:N)'; cn=(2*n+l)./(2*n.*(n+l»; 

Ba=B*a; [Hn,DHn]=shan3(Ba,N); 

Hn = Hn.'; DHn = DHn.'; 

Dn=DHn/Ba; an=(cn.*In)./(j*eta*Dn); 

Ets=j*eta*Pnl *(an.*DHn)/Ba; Ett=Eti+Ets; % ideally Ett --> 0 
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Hps=Pn 1 *(an. *Hn)/Ba; Hpt=Hpi+Hps; % Jtheta=H_pt 

Jtheta = -Hpt; 

Br=B*r; [Hn,DHn]=shan3(Br,N); 

Hn = Hn.'; DHn = DHn.'; 

Ets=j*eta*Pnl *(an.*DHn)/Br; Hps=Pnl*(an.*Hn)/Br; 
numring = Nth; % number of rings 
inc = input('Enter number of face segments: '); 
ppoints = inputCEnter phipoints - wavelength factor: '); 

% Field Point Inputs 

% Input the Field Point in terms of theta, 
radius = r; % radius of field point 
numFP = length(Hps); 
thetaFP = linspace(0,pi,numFP); 
thetaFPd= 180*thetaFP/pi; 
rho = radius. *sin(thetaFP); 
z = radius. *cos(thetaFP); 

% Jt is the tangential surface current at the nodes. 

% The nodes are defined by rhop & zp. 

Jt = Jtheta; 

% Determine the thetas for the sphere 

dtheta = pi/numring; % gives the delta theta 

theta = 0:dtheta:pi; % gives numring +1 values of theta. 

thetad = theta. * 1 80/pi; % theta in degrees 

% Determine (rhoprime,zprime) pairs on the sphere 
rhop = a*sin(theta); 
zp = a*cos(theta); 

% Find the length "s" of each face. 

s = sqrt(rhop(2) A 2 + (zp(l)-zp(2)) A 2); 

% Find the angle associated with each face 

alpha = acos(rhop(2)/s) + theta(l .numring); % radians 
alphad = alpha* 1 80/pi; % degrees 

talpha = tan(alpha); 
alphal = ones(inc+l,l)*alpha; 
alpha = reshape(alphal,numring*(inc+l),l); 
dzf = (zp(l :numring)-zp(2:numring+l ))/inc; 
phipoints = round(pi.*a*ppoints/lambda)+l; 
dphi = pi./phipoints; 

phi =0:dphi:pi; 

cp = cos(phi); 

for n = l:numring; % loops through rings 
zf(:,n) = (zp(n) : -dzf(n):zp(n+ 1 ))'; 
rhof(:,n) = (linspace(rhop(n),rhop(n+l),inc+l))'; 

JT(:,n) = Jt(n) *(zf( 1 : inc+ 1 ,n)-zp(n+ 1 ))/(zp(n)-zp(n+ 1 ))+. . . 

Jt(n+ 1 )*(zf( 1 :inc+l ,n)-zp(n))/(zp(n+l )-zp(n)); 
end 
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% Reshape the matrices to column vectors 
zf = reshape(zf,numring*(inc+l),l); 
rhof = reshape(rhof,numring*(inc+l), 1); 

JT = reshape(JT,numring*(inc+l),l); 
for FP = 1 :numFP; % loops through field points 
R1 = 2*rho(FP)*rhof*cp; 

R2 = (rho(FP) A 2+rhof. A 2+(z(FP)-zf) A 2)*ones(l,length(cp)); 

R = sqrt(R2-Rl); 

[Fes, FIs] = fcsfls2(cp,dphi,R); 

Ha = rhof.*sin(alpha); 

Fib = (z(FP) - zf).*cos(alpha); 

He = (Ha-Hb).*Fcs; 

Hd = rho(FP).*sin(alpha).*Fls; 

H = (Hc-Hd).*JT; 

H = reshape(H,inc+l,numring); 

HI = a*((l/(4*pi))*(tr apz(H))) ; 

HI = dzf.*Hl; 

HH(l:numring,FP) = HI.'; 
end 

Hphi = sum(HH); 

%% RMS Error Calculations 

E = 100*sqrt((Hps - Hphi.')'*(Hps - Hphi.'))./sqrt(Hps'*Hps); 
figure 

pIot(thdeg,abs(Hphi),'ko',thdeg,abs(Hps),'k') 

title([Hps (exact: integrated: o): ',eval('case'),'; RMS Error= ',num2str(E),'%']) 

xlabel('Theta (degrees)'); 

ylabel('|Hscattered|'); 

figure 

plot(thdeg,real(Hphi),'k*',thdeg,imag(Hphi),'k+',thdeg,real(Hps), . . . 
'k',thdeg,imag(Hps),'ko'); 

title([Hexact(reaI: imag: o) & Hint(real: *, imag: +): ',eval('case')]) 

xlabel('Theta (degrees)') 
ylabel(Hps (Amps/meter)') 

%text(. 55,. 95, ['Frequency= ',num2str(fO),'Mhz'], 'color', . . . 

% 'g', 'units', 'normalized'); 

%text(.55,.90,['Wavelength=',num2str(lambda),'meters'],'coIor',... 

% 'g', 'units', 'normalized'); 

%text(.55,. 85, [Dipole Offset^ ',num2str(z01),'*lambda'],'color',... 

% 'g','units','normalized'); 

%text(.55,.80,['Radius of Sphere= ',num2str(al),'*lambda'], 'color',... 

% 'g','units', 'normalized'); 

%text(.55,. 75, ['Radius of Field Point= ',num2str(rl),'*lambda'], 'color',... 

% 'g','units','normalized'); 

%text(.55,.70,[ r Number Rings= ’,num2str(numring)], 'color',... 

% 'g','units','normalized'); 
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%text(.55,.65,[Tace Segments= ',num2str(inc)],'color',... 

% 'gVunits','normalized'); 

%text(.55,. 60, ['Phi-Wavelength Factor= ',num2str(ppoints)],'color',... 

% 'g','unitsVnormalized'); 

%text(.55,.55, ['Field Points= ',num2str(numFP)],'color',... 

% 'g', 'units', 'normalized'); 

%text(.55,.50,[Tield Points per wl= ',num2str(ppwl)], 'color',... 

% 'g', 'units', 'normalized'); 

%text(.55,. 45, ['Number modes= ',num2str(N)],'color',... 

% 'g', 'units', 'normalized'); 

%text(.55,.40,['RMS Error= ',num2str(E),' %'], 'color',... 

% 'r','units', 'normalized'); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Computes Hphi incident at the field point radius due to the dipole. 

% This result is added to the Hphi scattered computed using Dan's 
% integration program 
dth=pi/Nth; th=(0:dth:pi)'; thdeg=180*th/pi; 
sth=sin(th); cth=cos(th); 

% computing rd, sin(thd), and cos(thd) for sphere surface (r=a) 
rd=sqrt(r*r+zO*zO-2*r*zO*cth); sthd=r*sth./rd; cthd=(r*cth-zO)./rd; 

% computing field components in dipole centered coordinates 
Fth=Idl*exp(-j*B*rd)./(4*pi*rd); 

Hpi=Fth.*sthd.*(j*B + l./rd); 

Hphi = Hphi.'; % make Hphi a column vector 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% This program loads the transfer function Hs from femdat.mat 
% and created by FEM2.M and dots it with Hphi stored in 
% source.mat from Dan's integration program SOURCE.M. 

% 

% The result is the current on the surface of the sphere 
% arrived at through the back propagation of the scattered 
% H field. 

% 

% The result (Jback) is compared to the surface current used in 
% Dan's program (Jtheta, stored in jtheta.mat). 

% 

% The surface current is generated by DIPSPHR2.M 
% 

load femdat.mat 

% Add Hphi scattered from sphere and Hpi incident from dipole 
Hphitot = -(Hphi + Hpi); 

% Determine Surface current due to back-propagation of Hphi 
% Hhpi is UNFILTERED 
Jback = Hs*Hphitot; 

% Determine Surface current due to back-propagation of Hphi 
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% RMS Error Calculations 

E = 100*sqrt((Jtheta - Jback)'*(Jtheta - Jback))./sqrt(Jtheta'*Jtheta); 
figure 

plot(thdeg,abs(Jtheta),'k',thdeg,abs(Jback),'ko') 
title(['Jexact (-) & Jback (o): ',eval('case'),'; RMS Error= ',... 

num2str(E),'%']) 
xlabel('Theta (degrees)') 
ylabel('|Jtheta|') 

%legend('Jtheta - dipsphr2.m', 'Jback: no filtering - fem2.m') 
figure 

plot(thdeg,real(Jtheta),'k',thdeg,imag(Jtheta),'ko',thdeg,real(Jback) > ... 

'k*',thdeg,imag(Jback),'k+') 

title(['Jexact(real: imag: o) & Jback(real: *, imag: +): ',eval('case')]); 

xlabel('Theta (degrees)') 
ylabel('Jtheta (Amps/m*m)') 
yn = input('Save data from this run? (Y/N): ','s'); 
if ~isempty(yn) 
if yn = y | yn = 'Y' 

clear eta Idl B a zO dth th sth cth rd sthd cthd Fth Etd Erd cdth sdth ... 

Eti r Sn Pnl In n cn Ba Hn DHn Ets Br Dn Ett FIs FP Fes E 
clear H HI HH Ha Hb He Hd JT Jt R R1 R2 alpha alphal alphad an cp dphi... 

dtheta dzf rb fho rhof rhop rs s talpha yn z zb zf zp zs 
clear Hpt Nmax Nth f phi phipoints... 

radius rho theta thetaFP thetaFPd thetad 
save (eval('case')) 

%save fO lambda al zOl rl Nth N inc ppoints Hs Hps numring ppoints numFP... 

% ppwl name thdeg Hps Jback Hphi Hphitot 

end 

end 
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function [Fc_s,Fl_s] = fcsfls2(cp,dphi,R) 

% For use with HPHIERR2.M 
% 

% This function does the integration to find Fc(s) & Fl(s) 

% updated 4/19/79 at home Dan Wawrzyniak 

lambda = 1; 

beta = 2* pi/lambda; 

phase = exp(-i*beta.*R); 

same = (l+i*beta.*R)./R. A 3; 

FI = same.*phase; 

FIs = (2*dphi*(trapz(Fl. '))).'; 

Fc = ones(size(Fl_s))*cp.*same.*phase; 

Fc_s = (2*dphi*(trapz(Fc. '))).'; 
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% DipSpr2.M computes scattered fields due to an elemental z-axis 
% dipole positioned above the north pole of a metallic sphere. 

% Mod-2 by M.A. Morgan 3/14/97 
% shan3.m & legpol2.m loaded 5/2/97 
clear all; 

eta=120*pi; % ffee-space Z_0 

Idl=l; % dipole moment 
fD=300; % 300 MHz 

B=2*pi*fD/300; % wavenumber 

zO=input('Enter z-position of elemental dipole in meters: '); 
a=input('Enter metal sphere radius in meters: '); 

Nth=in put ('Enter number of theta segments: '); 
dth=pi/Nth; th=(0:dth:pi)'; thdeg=180*th/pi; 
sth=sin(th); cth=cos(th); 

% computing rd, sin(thd), and cos(thd) for sphere surface (r=a) 
rd=sqrt(a*a+z0*z0-2*a*z0*cth); sthd=a*sth./rd; cthd=(a*cth-zO)./rd; 
% computing field components in dipole centered coordinates 
Fth=Idl*exp(-j*B*rd)./(4*pi*rd); 

Hpi=Fth.*sthd.*(j*B + 1 -/rd); 

Etd=eta*Fth.*sthd.*(j*B + 1 ./rd - j./(B*rd.*rd)); 

Erd=2*eta*Fth. *cthd. *(1 ./rd - j./(B *rd. *rd)); 

% translating incident field spherical components 
cdth=cth.*cthd + sth.*sthd; sdth=sth.*cthd - cth.*sthd; 

Eti=Etd.*cdth - Erd.*sdth; % note that Hpi=Hpd 
% Computing spherical harmonic coefficients for scattered field 
% Ets expansion which minimizes the total sphere surface 
% tangential field = Eti+Ets 
N=input('Enter number of spherical harmonics: '); 

Sn=-dth*Eti.*sth; Pnl=legpol2(cth,N); In=(Pnl.')*Sn; 
n=(l:N)'; cn=(2*n+l)./(2*n.*(n+l)); 

Ba=B*a; [Hn,DHn]=shan3(Ba,N); 

Hn = Hn.'; DHn = DHn.'; 

Dn=DHn/Ba; an=(cn.*In)./(j*eta*Dn); 

Ets=j*eta*Pnl*(an.*DHn)/Ba; Ett=Eti+Ets; % ideally Ett — > 0 

Hps=Pnl*(an.*Hn)/Ba; Hpt=Hpi+Hps; % Jtheta=H_pt 

Jtheta = -Hpt; 

save dipsphr2 Jtheta thdeg 

% Plotting expansions on sphere surface 

elf reset; plot(thdeg,abs(Eti),thdeg,abs(Ett),'.'); 

xlabel('theta (deg)'); 

ylabel('|E_theta|: Incident (solid); Total (dots)'); 

title(['Dipole at z0-,num2str(z0),'; Sphere Radius-, num2str(a),... 

'; No. Modes-, int2str(N)]); 
figure(l); 

hcpy=input('Enter 1 for Hard Copy: '); 
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if hcpy — 1, print; end; 

elf reset; plot(thdeg,abs(Hpi),thdeg,abs(Hps),'.',thdeg,abs(Hpt),'~'); 
xlabel('theta (deg)'); 

ylabel('|H_phi|: Incident (solid); Scattered (dots); Total (dashed)'); 
title([T>ipole at z0=',num2str(z0),'; Sphere Radius-, num2str(a),... 

'; No. Modes- ,int2str(N)]); 
figure(l); 

hcpy=input('Enter 1 for Hard Copy: '); 
if hcpy = 1, print; end; 

% Now compute scattered fields at specified radius 
r=input('Enter radius in meters to field point: '); 

Br=B*r; (Hn,DHn]=shan3(Br,N); 

Hn = Hn.'; DHn = DHn.'; 

Ets=j*eta*Pnl *(an.*DHn)/Br; Hps=Pnl *(an.*Hn)/Br; 

% Plotting scattered field expansions 

elf reset; plot(thdeg,abs(Ets),thdeg,eta*abs(Hps),'.'); 

xlabel('theta (deg)'); 

ylabel('|E_theta| (solid); eta*|H_phi| (dots)'); 
title(['Scattered Fields at r=',num2str(r),' m; Dipole at zO-,... 
num2str(z0),'; Sphere Radius- ,num2str(a),'; No. Modes- ,int2str(N)]); 
figure(l); 

hcpy=input('Enter 1 for Hard Copy: '); 

if hcpy = 1, print; end; 

end; 
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% FEM2.M provides finite-element solution for user-specified 
% metallic body of revolution using VARELA algorithm. 

% By M.A. Morgan 3/20/97- 4/17/97 
clear all 

% Generating surface and boundary points then mesh topology 
sgen=input(Enter name of surface generation program: 
eval(['[rs zs rb zb]-,sgen,';']); 

[rho zee nd elnd]=mesh2(rs,zs,rb,zb); % Computing mesh 
% Computing Mesh Parameters 
N=length(rho); % Number of Nodes in Mesh 
L=length(elnd(:, 1)); % Number of Elements in Mesh 

Ns=length(nd); % No. of Surface or Boundary Nodes 
Nm=N-nd(l)-nd(Ns)-Ns+2; % Number of Nodes for Unknown Hphi 
disp(’ ’) 

disp(Tinite Element Mesh Parameters: ') 
disp([' Number of Nodes: ’,int2str(N)]) 

disp([' Number of Elements: ',int2str(L)]) 
disp([' Number of Unknowns: ’,int2str(Nm)]) 
dispC ') 

disp(Tress a Key to Display Mesh or <Ctl> C to Abort'); 
pause 

% Plotting Mesh Nodes and Elements 
%clf reset 

figure('PaperPosition',[1.5 0.5 5 3.75]) 
if Ns <=19, plot(rho,zee,'go'); 

else, plot(rho,zee,'k.'); end 
hold on 
for 1=1 :L 
nds=elnd(l,:); 

rl=[rho(nds); rho(nds(l))]; zl=[zee(nds); zee(nds(l))]; 
plot(rl,zl,'k'); 
end 

v=axis; v(l)=0; v(2)=v(4)-v(3); axis(v); axis square; 
title('2-D Axisymmetric Shape with mesh') 
xlabel('meters') 
ylabel('meters') 

text(2.7*v(2)/5,5*v(4)/6,'Mesh Parameters: ') 

text(3 *v(2)/5,3 *v(4)/4, ['Nodes: ’,int2str(N)]) 

text(3 *v(2)/5,2*v(4)/3, [’Elements: ’,int2str(L)]) 

text (3 *v(2)/5, 7 *v(4)/ 12, [’Unknowns: ’,int2str(Nm)]) 

yn=input(’Print Hard Copy ? (Y/N): ’,'s’); 

if ~isempty(yn), if yn = 'y' | yn = 'Y', print; end; end 

% Loading System Matrices By Indexing Through Elements 

disp('Sparse Matrix Loading Using Element Contributions Can Now Begin') 

f=input('Enter Frequency in MHz to Start Loading; <Ctl> C to Stop: '); 
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k2=(pi*f7l50) A 2; T=zeros(3,3); 

A=sparse(Nm,Nm); B=sparse(Nm,Ns-2); 

12=0; nf=cumsum(nd); ns=nf-nd+l; % node finish/start per row 
for i=2:Ns; 

11=12+1; 12=ll+nd(i)+nd(i-l)-3; % start/finish elements 
for 1=11:12 

nds=elnd(l,:); rp=rho(nds); zp=zee(nds); 

T=varela(k2,rp,zp); 

disp(['Loading Element ',int2str(l),' of ',int2str(L)]) 

% Loading sparse matrix contributions 
for m=l:3; 
ma=0; 

if nds(m) ~= nf(i- 1 ) & nds(m) ~=nfl(i), % m not on boundary 
if nds(m) > nf(i-l) & i < Ns, 
ma=nds(m)-ns(2)-i+3; end 
if nds(m) < ns(i) & i > 2, 
ma=nds(m)-ns(2)-i+4; end 
if ma >0, % m is a solution node 

for n=l:3; 
na=0; 

if nds(n) ~= nf(i-l) & nds(n) ~=nf(i), 
if nds(n) > nf(i-l) & i < Ns, 
na=nds(n)-ns(2)-i+3; end 
if nds(n) < ns(i) & i > 2, 
na=nds(n)-ns(2)-i+4; end 
if na > 0, 

A(ma,na)= A(ma,na)+T(m,n); end 
end % n not on boundary end 
if nds(n) == nf(i-l) & i > 2, 

B(ma,i-2)=B(ma,i-2)-T(m,n); end 
if nds(n) = nf(i) & i < Ns, 

B(ma,i-l)=B(ma,i-l)-T(m,n); end 
end % n-loop end 
end % ma > 0 end 
end % m not on boundary end 
end % m-loop end 
end % 1-loop end 
end % i-loop end 

dispCMatrices Loaded, Now Solving System ... Be Patient') 

% Solve sparse matrix system to construct transfer function 
% array relating internal Hphi to Ns-2 nodal boundary values 
H = A\B; % clear A B C D % freeing memory 
% Extract rows to form transfer array relating Ns nodal values 
% of Hphi (including two z-axis Hphi=0 BC's) to Ns nodal values 
% of Hphi on the PEC surface. Note that the PEC surface output 
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% column vector has Hphi(l)=Hphi(Ns)=0 since these correspond 
% to nodes on the z-axis. 

Hs=zeros(Ns,Ns); 

row=l; Hs(2,2:Ns-l) = H(row,:); 

for n=3:Ns-l; 

row=row+nd(n-l)-l; Hs(n,2:Ns-l) = H(row,:); 
end 

n=l:Ns; elf; plot(n,diag(Hs),'r') 

v=axis; axis([l Ns v(3) v(4)]) 

xlabel('row(n)'); ylabel('Real diag[Hs]') 

title(['Diag[Hs] for ',sgen,'; Ns-,int2str(Ns),'; N-,int2str(N),... 

Nm-,int2str(Nm),'; L-,int2str(L)]); figure(l) 
yn=inputCPrint Hard Copy ? (Y/N): 
if ~isempty(yn), if yn = 'y' | yn = 'Y, print; end; end 
% Saving Needed Parameters and Hs Transfer Function Array 
save femdat f rs zs rb zb Hs 
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% Program FEMCHK4.M compares exact and FEM2.M computed 
% nodal values for spherical harmonic single mode solutions 
% of an offset metal sphere. By M.A. Morgan 4/21/97 
% Mod-2a (4/30/97) uses MatLab Bessel and Hankel Functions in shan3.m 
% Mod-3 (5/2/97) uses Wronskian to simplify r=a "exact" solution 
% and computes and plots errors for applied modes 
% Mod-4 (5/3/97) uses real mode BC's Pn A l(cos(theta)) on r=b 
% which produces a real internal solution 
clear all; 

case = input('What case are you running? ','s'); 
load femdat; % loading f rs zs rb zb Hs 

k=pi*f/150; Ns=length(rs); 

d=(zs(l)+zs(Ns))/2 - (zb(l)+zb(Ns))/2; % computing offset 

a=zs(l)-d; b=zb(l); Ra=k*a; Rb=k*b; % reset sphere centers 
zbd=zb-d; r_b=sqrt(rb. A 2+zbd. A 2); b=zb(l); 
th=linspace(0,pi,Ns)'; c_a=cos(th); th_a=linspace(0, 1 80, Ns)'; 
c_b=zbd./r_b; th_b=l 80*acos(c_b)/pi; 

N=input('Enter Upper Spherical Harmonic Mode to Check: '); 
[hna,Dhna]=shan3 (Ra,N); [hnb,Dhnb]=shan3 (Rb,N); 

DJna=real(Dhna); Jnb=real(hnb); 

DNna=-imag(Dhna); Nnb=-imag(hnb); 

Dn=(Jnb.*DNna-DJna.*Nnb); An=ones(Ns,l)*DNna; Bn=ones(Ns,l)*DJna; 
[Hn,DHn]=shan3(k*r_b,N); Jn=real(Hn); Nn=-imag(Hn); 
pna=legpol2(c_a,N); pnb=legpol2(c_b,N); pct=zeros(N,l); 
Ha=b*pna./(a*ones(Ns,l)*Dn); % "exact" r=a 
Hb=b*pnb.*(An.*Jn-Bn.*Nn)./(r_b*Dn); % "exact" r=b 
Hc=Hs*Hb; % Computed solution on r=a 

for n=l:N; 

pct(n)=100*sqrt((Ha(:,n)-Hc(:,n))'*(Ha(:,n)-Hc(:,n)))/... 

sqrt(Ha(:,n)'*Ha(:,n)); 

end 

[tsegl,tseg2] = size(Hs); 
figure('PaperPosition',[1.5 0.5 5 3.75]) 
bar(l:N,pct,'k'), v=axis; v(l)=0; v(2)=N+l; v(3)=0; axis(v); 
xlabel(Mode Index "n"'); ylabel('Percent'); 

title(['Single Mode RMS Error for Centered PEC Sphere: ',eval('case')]) 
yn=input('Print Hard Copy ? (Y/N): ','s'); 
if ~isempty(yn), if yn = 'y' I yn = 'Y', print; end; end 
while 1, 

n=input(['Enter n <= ',int2str(N),... 

' for Spherical Harmonic (0 to stop): ']); 
if isempty(n), break; 
elseif n = 0, break; 

elseif n > N, disp('n exceeds N'); break; end 
figure('PaperPosition',[1.5 0.5 5 3.75]) 
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plot(th_b,Hb(:,n),’k') 

xlabel('Theta (deg)'); ylabelCHphi (A/m)’); 
title(['Single Mode Hphi at Outer Boundary: ',eval('case'),... 
n=',int2str(n)]) 

v=axis; axis([0 180 v(3) v(4)]); 

yn=inputCPrint Hard Copy ? (Y/N): 

if ~isempty(yn), if yn = 'y' | yn = 'Y', print; end; end 

if length(th_a) <38, 

figure(PaperPosition',[1.5 0.5 5 3.75]) 

plot(th_a,Ha(:,n),'k’,th_a,Hc(:,n),'k.') 

title([Exact (line), FEM (dots): ',eval(’case'),... 

'; n-,int2str(n),'; RMS Error=',num2str(pct(n)),'%']) 

else, 

figureCPaperPosition',[1.5 0.5 5 3.75]) 
plot(th_a,Ha(:,n),'k',th_a,Hc(:,n),'k.'); 
title([Exact (line), FEM (dots): ',eval('case'),... 

’; n=',int2str(n),'; RMS Error=',num2str(pct(n)),'%']) 
end 

xlabel('Theta (deg)'); ylabel('Hphi (A/m)'); 
v=axis; axis([0 180 v(3) v(4)]); 
yn=inputCPrint Hard Copy ? (Y/N): ','s'); 
if ~isempty(yn), if yn = 'y' | yn = 'Y', print; end; end 
end 
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% hspherel.m 

% Program Hsphere.M uses spherical harmonics to compute 
% the Hs array for a centered PEC sphere. Saves in same 
% format as FEM2.m and can be used in FEMCHKxx.m programs 
% By M.A. Morgan 5/2/97 
% 5/10/97 modified by Dan Wawrzyniak 
% Nmax = 2*k*b 
clear all 

fHnputCEnter frequency in MHz: '); k=pi*f/150; lambda = 3e2/f; 
a=input('Enter metal sphere radius "a" in meters: '); 
b=inputCEnter outer mesh radius "b" in meters: '); 

Nsl=input('Enter number of theta points per wavelength: '); 

Ns = fix(pi*b*Nsl /lambda); 
disp(['Total theta points = ',int2str(Ns)]) 

th=linspace(0,pi,Ns); thd=linspace(0,180,Ns); dthl2=pi/(12*(Ns-l)); 
Ra=k*a; Rb=k*b; cth=cos(th); sth=sin(th); Nmax=2*fix(Rb); 
disp(' '); disp(['Estimated Nmax=integer[2*k*b] is ',int2str(Nmax)]) 
N=input(TEnter Upper Mode Order to Use in Computing Hs: '); 
disp(' ') 

[Hna,DHna]=shan3 (Ra,N); [Hnb,DHnb]=shan3(Rb,N); 
Pnl=legpol2(cth,N); n=l:N; Cn=(2*n+l)./(2*n.*(n+l)); 
Qn=(b/a)./(real(DHna). *imag(Hnb)-real(Hnb). *imag(DHna)); 
n=l:N; A=zeros(N,Ns); Hs=zeros(Ns,Ns); 

% Numerical integration of up(th)*Pnl(cth)*sth 
for p=2:Ns-l; 

fh=Pnl(p-l:p+l,:); g=sth(p-l:p+l); 
Inp=dthl2*(fh(l,:)*g(l)+6*fh(2,:)*g(2)+fh(3,:)*g(3)+... 

fi 1 (l,:)*g(2)+fii(2,:)*g(l)+fh(2,:)*g(3)+fh(3,:)*g(2)); 

A(:,p)=(Cn.*Qn.*Inp).'; 

end 

Hs=Pnl*A; 

n=l :Ns; elf; plot(n,diag(Hs),'r') 
v=axis; axis([l Ns v(3) v(4)]) 
xlabel('row(n)'); ylabel('Real Hs Diagonal') 
title(['Diag[Hs] using Hsphere.m for Ns=',int2str(Ns),... 

Nmax=',int2str(N),]); 

figure(l) 

yn=input('Print Hard Copy ? (Y/N): ','s'); 
if ~isempty(yn), if yn = 'y* | yn = 'Y', print; end; end 
rs=a*sth'; zs=a*cth'; rb=b*sth'; zb=b*cth'; 

% Saving Needed Parameters and Hs Transfer Function Array 
save femdat f rs zs rb zb Hs 
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function [Hn,DHn] = shan3(X,N) 

% Riccati Spherical Hankel Function Hn A (2)(X) = Jn(X) -j*Yn(X) 

% and Derivative. Using MatLab cylindrical functions of order n+1/2. 
% Function returns: Hn(n)=Hn A (2)(x) and DHn(n)=dHn A (2)(x)/dx 
% for n=l to N. Input N is a positive integer. X is an array 
% of real or complex values [xl x2 ... xM], Output Hn(xm) and 
% DHn(xm) are M by N complex arrays. By M. A. Morgan. 

% Mod-3 (5/1/97) allows X to be either row or column array. 

[ml m2]=size(X); M=max(ml,m2); Hn=zeros(M,N); DHn=Hn; 
if ml=l, x=X*; else, x=X; end 
Nu=(l:N)+0.5; xn=x*ones(l,N); 

J0=sqrt(pi*x/2).*besselj(0.5,x); 

Jn=sqrt(pi*xn/2).*besselj(Nu,x); 

N0=sqrt(pi*x/2).*bessely(0.5,x); 

Nn=sqrt(pi*xn/2).*bessely(Nu,x), 

Hn=Jn-j*Nn; H0=J0-j*N0; DHn(:,l)=HO-Hn(:,l)./x; 
ifN > 1, DHn(:,2)=Hn(:,l)-2*Hn(:,2)./x; end 
if N > 2, 

for n=3:N; DHn(:,n)=Hn(:,n-l)-n*Hn(:,n)./x; end 
end 
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function Pnl = legpo!2(X,N) 

% Computing a matrix of associated Legendre polynomials 
% Pn A m(x) with m=l, for n=l, 2, .. N where N is a positive integer 
% and X is an M-element array of real values [xl x2 ... xM], 

% Computed real array Pnl(xm) has M rows and N columns. 

% Using upward recurrence formula from page 953 of Balanis - 

% Advanced EM Engineering, Wiley, 1989. Program by M.A. Morgan 

% Mod-2 (5/1/97) allows either row or column input X array 

[ml m2]=size(X); M=max(ml,m2); Pnl=zeros(M,N); 

if m 1=1, x=X'; else, x=X; end 

Pnl =zeros(M,N); Pn 1 ( : , 1 )=- sqrt( 1 -x. *x); 

ifN > 1, Pnl(:,2)=3*x.*Pnl(:,l); end 

if N > 2, 

for n=2:N-l; 

nl=l/n; Pnl(:,n+l)=(2+nl)*x.*Pnl(:,n)-(l+nl)*Pnl(:,n-l); 

end 

end 
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function [rs,zs,rb,zb] = osphere 

% Computing meridian surface coordinates for offset sphere 
% (rs,zs) with a=radius and d=z-axis offset. Outer 
% mesh boundary is centered b=radius sphere (rb,zb). 

% Returning coordinates in Ns x 4 array. 

% By M. A. Morgan 3/24/97 

disp('Offset Sphere Surface and Mesh Boundary Program'); 
a=input('Enter sphere radius (meters): '); 
d=inputCEnter sphere offset (meters): '); 
b=inputCEnter mesh radius (meters): '); 

Ns=inputCEnter number of surface points: '); 
tha=linspace(0,pi,Ns)'; sta=sin(tha); cta=cos(tha); 
rs=a*sta; zs=a*cta+d; 

thb=tha-asin(d*sin(pi-tha)/b); stb=sin(thb); ctb=cos(thb); 
rb=b*stb; zb=b*ctb; 

elf reset; plo^rs.zs.rs^s/.g^rb.zb.rb^b/.g') 
v=axis; v(l)=0; v(2)=v(4)-v(3); axis(v); axis square; 
xlabel('Press a Key to Display Mesh or <Ctrl> C to Abort') 
figure(l) 
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function [rs,zs,rb,zb] = cone 

% Computing meridian surface coordinates for centered cone 
% (rs,zs) with h=height and b=base radius . Outer 
% mesh boundary is centered a=radius sphere (rb,zb). 

% Returning coordinates in Ns x 4 array. 

% By M.A. Morgan 3/27/97 

disp('ConicaI Surface and Spherical Mesh Boundary Program'); 
h=inputCEnter cone height (meters): ’); 
b=inputCEnter cone base radius (meters): '); 
a=inputCEnter outer mesh spherical radius (meters): '); 
Ns=inputCEnter number of surface points: '); 
rs=zeros(Ns,l); zs=rs; th=rs; 

S=b+sqrt(b*b+h*h); % total surface length on cone 

Nb=fix(Ns*b/S)+l; ifNb>Ns-2, Nb=Ns-2; end; db=b/Nb; Nz=Ns-Nb; 

rs(l :Nz)=linspace(0,b,Nz); zs(l :Nz)=linspace(h/2,-h/2,Nz); 

rs(Nz+l :Ns)=linspace(b-db,0,Nb); zs(Nz+l :Ns)=(-h/2)*ones(Nb, 1 ); 

% distort boundary node positions to improve mesh near comer 

a2=atan(h/b)/2; gm=a2+pi/2; c=(h/2)-b*tan(a2); 

dl=asin(c * sin(gm)/ a); thc=gm+dl; dtb=(pi-thc)/Nb; 

th(l :Nz)=linspace(0,thc,Nz); th(Nz+l :Ns)=linspace(thc+dtb,pi,Nb); 

st=sin(th); ct=cos(th); rb=a*st; zb=a*ct; 

elf reset; plot(rs,zs,rs,zs,'r.',rb,zb,rb,zb,'g.') 

v=axis; v(l)=0; v(2)=v(4)-v(3); axis(v); axis square; figure(l) 

xlabel('Press a Key to Display Mesh or <Ctrl> C to Abort') 
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